Variational method for finding periodic orbits in a general flow 
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A variational principle for determining unstable periodic orbits of flows as well as unstable spatio- 
temporally periodic solutions of extended systems is proposed and implemented. An initial loop 
approximating a periodic solution is evolved in the space of loops toward a true periodic solution by 
a minimization of local errors along the loop. The "Newton descent" partial differential equation 
that governs this evolution is an infinitesimal step version of the damped Newton-Raphson iteration. 
The feasibility of the method is demonstrated by its application to the Henon-Heiles system, the 
circular restricted three-body problem, and the Kuramoto-Sivashinsky system in a weakly turbulent 
regime. 
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I. INTRODUCTION 



The periodic orbit theory of classical and quantum 
chaos 0, 01 is one of the major advances in the study 
of long-time behavior of chaotic dynamical systems. The 
theory expresses all long time averages over chaotic dy- 
namics in terms of cycle expansions 0, sums over pe- 
riodic orbits (cycles) ordered hierarchically according to 
the orbit length, stability, or action. If the symbolic dy- 
namics is known, and the flow is hyperbolic, longer cycles 
are shadowed by the shorter ones, and cycle expansions 
converge exponentially or even super-exponentially with 
the cycle length 

A variety of methods for determining all periodic orbits 
up to a given length have been devised and successfully 
implemented for low-dimensional systems 0, IE S B B - 
For more complex dynamics, such as turbulent flows [lOj, 
nonlinear waves [13, or quantum fields jTHIl with high 
(or infinite) dimensional phase spaces and complicated 
dynamical behavior, many of the existing methods be- 
come unfeasible in practice. In the most computationally 
demanding calculation carried out so far, Kawahara and 
Kida 01 have found two periodic solutions in a 15,422- 
dimensional discretization of a turbulent plane Couette 
flow. The topology of high-dimensional flows is hard to 
visualize, and even with a decent starting guess for the 
shape and location of a periodic orbit, methods like the 
Newton-Raphson method are likely to fail. In ref. |0 
we have argued that variational, cost-function minimiza- 
tion methods offer a robust alternative. Here we derive, 
implement and discuss in detail one such new variational 
method for finding periodic orbits in general flows, and 
specifically high-dimensional flows. 

In essence, any numerical algorithm for finding peri- 
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odic orbits is based on devising a new dynamical system 
which possesses the desired orbit as an attracting fixed 
point with a sizable basin of attraction. Beyond that, 
there is much freedom in constructing such system. 

For example, the multipoint shooting method elimi- 
nates the long-time exponential instability of unstable or- 
bits by splitting an orbit into a number of short segments, 
each with a controllable expansion rate. The multiple 
shooting combined with the Newton-Raphson method is 
an efficient tool for locating periodic orbits of maps [T6l |. 
A search for periodic orbits of a continuous time fiow can 
be reduced to a multiple shooting search for periodic or- 
bits of a set of maps by constructing a set of phase space 
Poincare sections such that an orbit leaving one section 
reaches the next one in a qualitatively predictable man- 
ner, without traversing other sections along the way. In 
turbulent, high-dimensional flows such sequences of sec- 
tions are hard to come by. One solution might be a large 
set of Poincare sections, with the intervening flight seg- 
ments short and controllable. 

Here we follow a different strategy, and discard 
Poincare sections altogether; we replace maps between 
spatially flxed Poincare sections, by maps induced by 
discretizing the time evolution into small time steps. For 
sufficiently small time steps such maps are small defor- 
mations of identity. We distribute many points along a 
smooth loop L, our initial guess of a cycle location and 
its topological layout. If both the time steps and the loop 
deformations are taken infinitesimal, a partial differential 
equation governs the "Newton descent" , a fictitious time 
flow of a trial loop L into a genuine cycle p, with expo- 
nential convergence in the flctitious time variable. We 
then use methods developed for solving PDEs to get the 
solution. Stated succinctly, the idea of our method is to 
make an informed rough guess of what the desired cycle 
looks like globally, and then use a variational method to 
drive the initial guess towards the exact solution. For ro- 
bustness, we replace the guess of a single orbit point by a 
guess of an entire orbit. For numerical safety we replace 
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the Newton-Raphson iteration by the "Newton descent" , 
a differential flow that minimizes a cost function com- 
puted as deviation of the approximate flow from the true 
flow along a smooth loop approximation to a cycle. 

In sect. |n]we derive the partial differential equation 
which governs the evolution of an initial guess loop to- 
ward a cycle and the corresponding cost function. An 
extension of the method to Hamiltonian systems and sys- 
tems with higher time derivatives is presented in scct. lIIII 
Simplifications due to symmetries and details of our nu- 
merical implementation of the method are discussed in 
sect. IIVI In sect. [3 we test the method on the Henon- 
Heiles system, the restricted three body problem, and 
a weakly turbulent Kuramoto-Sivashinsky system. We 
summarize our results and discuss possible improvements 
of the method in sect. I VII 



II. THE NEWTON DESCENT METHOD IN 
LOOP SPACE 



A. A variational equation for the loop evolution 

A periodic orbit is a solution {x,T), x e W^, T e R of 
the periodic orbit condition 



T> 



(1) 



for a given flow or discrete time mapping x h- > 

Our goal is to determine periodic orbits of flows deflned 

by flrst order ODEs 



dx 



v{x) , 



X eM c 



in many (even inflnitely many) dimensions d. Here A4 is 
the phase space (or state space) in which evolution takes 
place, TM is the tangent bundle and the vector field 
v{x) is assumed to be smooth (sufficiently differentiablc) 
almost everywhere. 

We make our initial guess at the shape and the location 
of a cycle p by drawing a loop L, a smooth, diffcrcntiable 
closed curve i(s) € L C Ai, where s is a loop parameter. 
As the loop is periodic, we find it convenient to restrict 
s to [0, 27r], with the periodic condition x{s) ~ x{s + 2tt). 
Assume that L is close to the true cycle p, pick N pairs 
of nearby points along the loop and along the cycle 

Xn = i(s„) , < Si < . . . < stv < 27r , 

Xn = X{tn) , < ti < . . . < tAT < Tp , (3) 

and denote by 5xn the deviation of a point Xn on the 
periodic orbit p from the nearby point Xn, 

The deviations 5x are assumed small, vanishing as L ap- 
proaches p. 



The orientation of the s-velocity vector tangent to the 
loop L 



v{x) 



dx 

ds 



is intrinsic to the loop, but its magnitude depends on the 
(still to be specified) parametrization s of the loop. 

At each loop point 5:„ G i we thus have two vectors, 
the loop tangent w„ = v{xn) and the flow velocity Vn = 
v{xn)- Our goal is to deform L until the directions of w„ 
and Vn coincide for all n = 1, . . . , TV, N oo, that is 
L ~ p. To match their magnitude, we introduce a local 
time scaling factor 



\{Sn) = Atn/ASn , 



(4) 



where As„ = Sn+i ~ Sm n — 1, . . . , iV — 1 , Asn — 2tt — 
{sn ~ Si), and likewise for At„. The scaling factor A(s„) 
ensures that the loop increment As„ is proportional to 
its counterpart At„ -I- Stn on the cycle when the loop L 
is close to the cycle p, with Stn ^ as L p. 

Let x{t) = /*(x) be the state of the system at time 
t obtained by integrating |(2Jl, and J{x,t) — dx{t)/dx{Q) 
be the corresponding Jacobian matrix obtained by inte- 
grating 



dt-"^-^' ~ dx. 



with J{x, 0) = 1 . (5) 



Since the point Xn = Xn + <5i„ is on the cycle. 



(x, v) e TM (2) Linearization 



f^x) ^x + v{x)5t , f{x + 6x) « x{t) + J(x, t)dx , 

of © about the loop point Xn and the time interval At„ 
to the next cycle point leads to the multiple shooting 
Newton-Raphson equation, for any step size Af„: 

SXn+l - J{Xn, Atn)SXn - Vn+lStn = /^*"(iri) - Xn+1 ■ 

(7) 

Provided that the initial guess is sufficiently good, the 
Newton-Raphson iteration of iQ generates a sequence of 
loops L with a decreasing cost function |0| 



N 



N 



XN+1 = Xi 



(8) 

The prefactor N/ (27r)^ makes the definition of consis- 
tent with (|13|l in the iV oo limit. If the flow is locally 
strongly unstable, the neighborhood in which the lin- 
earization is valid could be so small that the full Newton 
step would overshoot, rendering bigger rather than 
smaller. In this case the step-reduced, damped Newton 
method is needed. As proved in ref. , under conditions 
satisfied here, decreases monotonically if appropriate 
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L(co)=p 




FIG. 1: (a) An annulus L{t) swept by the Newton descent 
flow dx/dr, connecting smoothly the initial loop 1/(0) to the 
periodic orbit p = L{oo). (b) In general the loop velocity field 
v{x) does not coincide with A«(i); for a periodic orbit p, it 
does so at every x £ p. 



step size is taken. If infinitesimal steps are taken, de- 
crease of is ensured. Wc parametrize such continuous 
deformations of the loop by a fictitious time r. 

We fix Asn and proceed by St each step of the itera- 
tion, that is, multiply the right hand side of Q by St. 
According to the change of Atn with respect to r is 

9A / 



f:p{Sn,T)STAs„. As SXn = ^x{Sn,T)S7 



equal to Str, 

dividing both sides of Q) by St yields 

dXn+l 



dT 



J(Xn,Atn)^^ 

dT 



u 

fAt 



9^, N A 

l-^{Sn,T)As„ 



/^*"(i„)-5n+l. (9) 



In the N ^ oo limit, the stepsizes As„,At„ = 
0{1/N) 0, and we have 

J(Xn,Atn) « 1 -I- A(Xn)Atn , ~ in + V^At^ . 

Substituting into lO and using the scaling relation 
we obtain 



d^x 
dsdT 



dx d\ 



dT 



Xv 



(10) 



This PDE, which describes the evolution of a loop L{t) 
toward a periodic orbit p, is the central result of this 
paper. The family of loops so generated is parametrized 



by a; = x{s,t) G L{t), where s denotes the position 
along the loop, and the fictitious time r parametrizes 
the deformation of the loop, see Fig. ^a). We refer to 
this infinitesimal step version of the damped Newton- 
Raphson method as the "Newton descent" . 

The important feature of this equation is that a de- 
creasing cost functional exists. Rewriting H1U|I as 



d 

— (w — Xv) = — (v — Xv) 

OT 



(11) 



we have 



w - Au = e"^(u - Aw)|r=o , (12) 
so the fictitious time r flow decreases the cost functional 



1 

2^ 



ds (v(x) 



Xv{x)f 



(13) 



L(r) 



monotonically as the loop evolves toward the cycle. 

At each iteration step the differences of the loop tan- 
gent velocities and the dynamical flow velocities are re- 
duced by the Newton descent. As t oo, the fictitious 
time fiow alignes the loop tangent v with the dynami- 
cal fiow vector v = Xv, and the loop x{s,t) g L(t), see 
Fig.^b), converges to a genuine periodic orbit p = L(oo) 
of the dynamical flow x = v{x). Once the cycle p is 
reached, by A(s,oo) — ^{x{s,oo)), and the cycle 
period is given by 



2ir 



A(i(s, oo))ds . 



Of course, as at this stage we have already idcntiflcd the 
cycle, we may pick instead an initial point on p and cal- 
culate the period by a direct integration of the dynamical 
equations 



B. Marginal directions and accumulation of loop 
points 

Numerically, two perils lurk in a direct implementation 
of the Newton descent (jlOII . 

First, when a cycle is reached, it remains a cycle under 
a cyclic permutation of the representative points, so on 
the cycle the operator 



A 



d_ 



XA 



has a marginal eigenvector v{x{s)) with eigenvalue 0. If A 
is flxed, as the loop approaches the cycle, H10() approaches 
its limit 

Ap = 0. 

OT 

Therefore, on the cycle, the operator A~^ becomes sin- 
gular and the numerical woes arise. 
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The second potential peril hides in the freedom of 
choosing the loop (re-)parametrization. Since s is re- 
lated to the time t by the yet unspecified factor A(s, r), 
uneven distributions of the sampling points over the loop 
L could arise, with the numerical discretization points .t„ 
clumping densely along some segments of L and leaving 
big gaps elsewhere, thus degrading the numerical smooth- 
ness of the loop. 

We remedy these difficulties by imposing constraints 
on H10() . In our calculation for Kuramoto-Sivashinsky 
system of sect. O the first difficulty is dealt with by in- 
troducing one Poincare section, for example, by fixing 
one coordinate of one of the sampling points, Xi(si, r) = 
const. This breaks the translational invariance along the 
cycle. Other types of constraints might be better suited 
to a specific problem at hand. For example, we can 
demand that the average displacement of the sampling 
points along the loop vanishes, thus avoiding a spiraling 
descent towards the desired cycle. 

We deal with the second potential difficulty by choos- 
ing a particularly simple loop parametrization. So far, 
the parametrization s is arbitrary and there is much 
freedom in choosing the best one for our purposes. We 
pick s— and r— independent constant scaling A(s, r) ~ A. 
With uniform grid size As„ = As and fixed A, the loop 
parameter s = t/X is proportional to time t, and the dis- 
cretization distributes the sampling points along the 
loop evenly in time. As the loop approaches a cycle, ^ 
is numerically obtainable from (|10|l . and on the cycle the 
period is given by Tp = 2ttX. 

Even though this paper focuses on searches for pe- 
riodic orbits, the Newton descent is a general method. 
With appropriate modifications of boundary conditions 
and scaling of time, (|10|l can be adapted to determination 
of homoclinic or heteroclinic orbits between equilibrium 
points or periodic orbits of a flow, or more general bound- 
ary value problems. Applied to 2-point boundary value 
problems, Newton descent is similar to the quasilineariza- 
tion ^3 but has the advantage that the free parameter 
A(s,t) is available for adjusting scales in the problem 
and that searches can be restricted to phase space sub- 
manifolds of interest. A simple example of a restriction 
to a submanifold are searches for cycles of a given en- 
ergy, constrained to the H{q,p) = E energy shell in the 
phase space of a Hamiltonian system. Furthermore, as 
we shall show now, the symplectic structure of Hamil- 
ton's equations greatly reduces the dimensionality of the 
submanifold that we need to consider. 



III. EXTENSIONS OF NEWTON DESCENT 

In classical mechanics particle trajectories are also so- 
lutions of a variational principle, the Hamilton's varia- 
tional principle. For example, one can determine a pe- 
riodic orbit of a billiard by wrapping around a rubber 
band of a roughly correct topology, and then moving 
the points along the billiard walls until the length (that 



is, the action) of the rubber band is extremal (maximal 
or minimal under infinitesimal changes of the boundary 
points). Note that the extremization of action requires 
only D configuration coordinate variations, not the full 
2£'-dimensional phase space variations. 

Can we exploit this property of the Newtonian me- 
chanics to reduce the dimenionality of our variational 
calculations? The answer is yes, and easiest to under- 
stand in terms of the Hamilton's variational principle 
which states that classical trajectories are extrema of the 
Hamilton's principal function (or, for fixed energy E, the 
action S = R + Et) 

R{qi,ti;qo,to) = [ dt £.{q{t) , q{t) , t) , 

Jto 

where C{q, q, t) is the Lagrangian. Given a loop L(t) we 
can compute not only the tangent "velocity" vector v, 
but also the local loop curvature or "acceleration" vector 

and indeed, as many s derivatives as needed. Matching 
the dynamical acceleration a(i) (assumed to be functions 
of X and v{x)) with the loop "acceleration" a{x) results 
in a new cost function and the corresponding PDE Hll|) 
for the evolution of the loop 

— (5 — X^a) ~ —(a — X^a) . 

OT 

We use A^ instead of A in order to keep the notation 
consistent with that is t = As. Expressed in terms 
of the loop variables x{s), the above equation becomes 

d^x ^da d^x ^2 / da dx \dX 

d'^sdr dv dsdr dx dr \dv ds J dr 

= X^a-a, (14) 

where v = Although (|14() looks more complicated 

than 11(J|) . in numerical fictitious time integrations, we 
are rewarded by having to keep only half of the phase 
space variables. 

More generally, if a differential equation has the form: 

x(") = /(x,x«,..- ,x('"-i)), (15) 

where x'^'^^ = , k = 1, - ■ ■ ,m and x e R**, the same 
technique can be used to match the highest derivatives 
A™a;(™) and i^™), 

OT 

with i*^™^ = -^-^x{s) calculated directly from i(s) on the 
loop by differentiation. In loop variables x{s) we have, 

d-'^'i ^„ y _dj_ ^y«-i2,(rn) 
QgmQj. dx^^^ dr X'^ds'^ dr 

= A™a;(™) - , (16) 
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where x = a;*^"^ and a;^'^^ = ;S^r^, k = I,-- - ,to — 1 
are assumed. Conventionally, 1)15(1 is converted to a sys- 
tem of md first order differential equations, whose dis- 
cretized derivative (see (|f 7|) below) are banded matrices 
with band width of bmd. Using Hlfci|l . we only need d 
equations for the same accuracy and the corresponding 
band width is (m + 4)d. The computing load has been 
greatly reduced, the more so the larger m is. Neverthe- 
less, choice of a good initial loop guess and visualization 
of the dynamics are always aided by a plot of the orbit 
in the full md-dimensional phase space, where loops can- 
not self-intersect and topological features of the flow is 
exhibited more clearly. 



IV. IMPLEMENTATION OF NEWTON 
DESCENT 

As the loop points satisfy a periodic boundary con- 
dition, it is natural to employ truncated discrete Fast 
Fourier Transforms (FFT) in numerical integrations of 
(|10|l . Since we are interested only in the final, stationary 
cycle p, the accuracy of the fictitious time integration is 
not crucial; all we have to ensure is the smoothness of the 
loop throughout the integration. The Euler integration 
with fairly large time steps 5t suffices. The computation- 
ally most onerous step in implementation of the Newton 
descent is the inversion of large matrix A in . When 
the dimension of the dynamical phase space of ^ is high, 
the inversion of A needed to get takes most of the inte- 
gration time, making the evolution extremely slow. This 
problem is partially solved if the finite difference meth- 
ods are used. The large matrix A then becomes sparse 
and the inversion can be done far more quickly. 



A. Numerical implementation 

In a discretization of a loop, numerical stability re- 
quires accurate discretization of loop derivatives such as 



dx 
ds 



{Di)r 



x—x[sji ) 



In our numerical work we use the four-point approxima- 
tion 13, 



/ 8 -1 
-8 8-1 
1-808 



D = 



12h 



1 -8\ 
1 



1 -8 8 -1 
1 1-808 
V 8 -1 1-80/ 

(17) 

where h = 2tt/N. Here, each entry represents a [dxd] 
matrix, 8 — > 81, etc., with blank spaces are filled with 
zeros. The two \2d x 2d\ matrices 



Ml = 



1 -81 
1 



Mo = 



-1 
81 -1 



located at the top- right and bottom-left corners take care 
of the periodic boundary condition. 

The discretized version of 1)1 n|l with a fictitious time 
Euler step St is 



A V 
d 



= 5i 



\v — V 




(18) 



where 



A = D + di'Ag[Ai,A2,--- ,An], 
with An = A{x{sn)) defined in Q, and 

V = {vi,V2r ■ ■ ,vn) , with u„ = u(x(s„)) , 

V = (ui,-U2, ■ ■ ■ ,'UAr) , with u„ = u(i(s„)) , 

are the two vector fields that we want to match every- 
where along the loop, a is an Nd dimensional row vector 
which imposes the constraint on the coordinate varia- 
tions Sx = ((5ii, (5x2, ■ ■ ■ ,Sxn). The discretized New- 
ton descent (I18|) is an infinitesimal time step variant of 
the multipoint (Poincare section) shooting equation for 
flows In order to solve for the deformation of the 

loop coordinates and period, Sx and SX, we need to in- 
vert the [{Nd+l)x{Nd+ 1)] matrix on the left hand 
side of CHJ. 

In our numerical work, this matrix is inverted using 
the banded LU decomposition on the embedded band- 
diagonal matrix, and the Woodbury formula on the 
cyclic, border terms. The LU decomposition takes most 
of the computational time and considerably slows down 
the fictitious time integration. We speed up the inte- 
gration by a new inversion scheme which relies on the 
smoothness of the flow in the loop space. It goes as fol- 
lows. Once we have the LU decomposition at one step, 
we use it to approximately invert the matrix in the next 
step, with accurate inversion achieved by the iterative 
approximate inversions pl] |. In our applications we flnd 
that a single LU decomposition can be used for many St 
evolution steps. The further we go, the more iterations at 
each step are needed to implement the inversion. After 
the number of such iterations exceeds some given fixed 
maximum number, we perform another LU decomposi- 
tion and proceed as before. The number of integration 
steps following one decomposition is an indication of the 
smoothness of the evolution, and we adjust accordingly 
the integration step size St: the greater the number, the 
bigger the step size. As the loop approaches a cycle, the 
evolution becomes so smooth that the step size can be 
brought all the way up to St = 1, the full undamped 
Ncwton-Raphson iteration step. In practice, one can 
start with a small but reasonable number of points, in 
order to get a coarse solution of relatively low accuracy. 
After achieving that, the refined guess loop can be con- 
structed by interpolating more points, and proceed with 
for a more accurate calculation in which St can be set as 
large as the full Newton step St — I, recovering the rapid 
quadratic convergence of the Ncwton-Raphson method. 
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It is essential that the smoothness of the loop is main- 
tained throughout the calculation. We monitor the 
smoothness by checking the Fourier spectrum of x{-,t). 
An unstable difference scheme for loop derivatives might 
lead to unbounded sawtooth oscillations l2a|. A heuristic 
local linear stability analysis (described in [23) indicates 
that our scheme is stable, and that the high frequency 
components do not generate instabilities. 



B. Initialization and convergence 

As in any other method, a qualitative understanding of 
the dynamics is a prerequisite to successful cycle searches. 
We start by numerical integration with the dynamical 
system Numerical experiments reveal regions where 
a trajectory spends most of its life, giving us the first 
hunch as to how to initialize a loop. We take the FFT 
of some nearly recurred orbit segment and keep only the 
lowest frequency components. The inverse Fourier trans- 
form back to the phase space yields a smooth loop that we 
use as our initial guess. Since any generic orbit segment is 
not closed and might exhibit large gaps, the Gibbs phe- 
nomenon can take the initial loop so constructed quite 
far away from the region of interest. We deal with this 
problem by manually deforming the orbit segment into 
a closed loop before performing the FFT. Searching for 
longer cycles with multiple circuits requires more delicate 
initial conditions. The hope is that a few short cycles can 
help us establish an approximate symbolic dynamics, and 
guess for longer cycles can be constructed by cutting and 
glueing the short, known ones. For low dimensional sys- 
tems, such methods yield quite good systematic initial 
guesses for longer cycles |2J|. 

An alternative way to initialize the search is by utiliz- 
ing adiabatic deformations of dynamics, or the homotopy 
evolution [2^ . If the dynamical system ^ depends on a 
parameter /i, short cycles might survive as /i varies pass- 
ing through a family of dynamical systems, giving in the 
process birth to new cycles through sequences of bifur- 
cations. Most short unstable cycles vary little for small 
changes of fi. So, a cycle existing for parameter value fii 
can be chosen as the initial trial loop for a nearby cycle 
surviving a small change /ii — > /X2. In practice, one or 
two iterations often suffice to find the new cycle. 

A good choice of the initial loop significantly expe- 
dites the computation, but there are more reasons why 
good initial loops are crucial. First of all, if we break the 
translational invariance by imposing a constraint such as 
ii{si,T) = c, we have to make sure that both the ini- 
tial loop and the desired cycle intersects this Poincare 
plane. Hence, the initial loop cannot be wildly different 
from the desired cycle. Second, in view of (|12(l . the loop 
always evolves towards a local minimum of the cost func- 
tional 113|l . with discretization points moving along the 
v — Xv fixed direction, determined by the initial condition. 
If the local minimum corresponds to a zero of the cost 
functional, we obtain a true cycle of the dynamical flow 



(0). However, if the value of cost functional is not equal 
to zero at the minimum while the gradient is zero, l|18|) 
yields a singular matrix A. In such cases the search has 
to be abandoned and restarted with a new initial loop 
guess. In the periodic orbit searches of sect. IVl starting 
with blind initial guesses (guesses unaided by a symbolic 
dynamics partition) , such local minima were encountered 
in about 30% of cases. 



C. Symmetry considerations 

The system under consideration often possesses certain 
symmetries. If this is the case, the symmetry should be 
both be feared for possible marginal eigen-directions, and 
be embraced as a guide to possible simplifications of the 
numerical calculation. 

If the dynamical system equations (O are invariant 
under a discrete symmetry, the concept of fundamental 
domain 5, 2GJ can be utilized to reduce the length of the 
initial loop when searching for a cycle of a given symme- 
try. In this case, we need discretize only an irreducible 
segment of the loop, decreasing significantly the dimen- 
sionality of the loop representation. Other parts of the 
loop are replicated by symmetry operations, with the full 
loop tiled by copies of the fundamental domain segment. 
The boundary conditions are not periodic any longer, but 
all that we need to do is modify the cyclic terms. Instead 
of using Ml and A/2 in H17() . we use MiQ and M2Q~^, 
where Q is the relevant symmetry operation that maps 
the fundamental segment to the neighbor that precedes 
it. In this way, a fraction of the points represent the cycle 
with the same accuracy, speeding up the search consid- 
erably. 

If a continuous symmetry is present, it may compli- 
cate the situation at first glance but becomes something 
that we can take advantage of after careful checking. For 
example, for a Hamiltonian system unstable cycles may 
form continuous families j27i Esj , with one or more mem- 
bers of a family belonging to a given constant energy sur- 
face. In order to cope with the marginal eigendirection 
associated with such continous family, we search for a 
cycle on a particular energy surface by replacing the last 
row of equation (jf 8(1 by an energy shell constraint [lfil |. 
We put one point of the loop, say X2, on the constant 
energy surface H{x) = E, and impose the constraint 
\jH{x2) ■ 5x2 = 0, so as to keep X2 on the surface for 
all T. The integration of H1U|) then automatically brings 
all other loop points to the same energy surface. Alter- 
natively, we can look for a cycle of given fixed period T 
by fixing A and dropping the constraint in the bottom 
line of 1(18(1 . These two approaches are conjugate to each 
other, both needed in applications. In most cases, they 
are equivalent. One exception is the harmonic oscillator 
for which the oscillations have identical period but dif- 
ferent energy. Note that in both cases the translational 
invariance is restored, as we have discarded the Poincare 
section condition of sect. IIIBl As explained in |^, this 
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causes no trouble in numerical calculations. 



V. APPLICATIONS 

We have checked that the iteration of (|18|) yields 
quickly and robustly the short unstable cycles for stan- 
dard models of low-dimensional dissipative flows such as 
the Rossler system A more daunting challenge are 
searches for cycles in Hamiltonian flows, and searches for 
spatio-temporally periodic solutions of PDEs. In all nu- 
merical examples that follow, the convergence condition 
is < 10-^. 



A. Henon-Heiles system and restricted three-body 
problem 

First, we test the Hamiltonian version of the Newton 
descent derived in Sect. IIIII by applying the method to 
two Hamiltonian systems, both with two degrees of free- 
dom. In both cases, our initial loop guesses are rather 
arbitrary combinations of trigonometric functions. Nev- 
ertheless, the observed convergence is fast. 

The Henon-Heiles system |33 is a standard model in 
celestial mechanics, described by the Hamiltonian 

H=^ipl+pl + x' + y^)+x^y-^. (19) 

It has a time reversal symmetry and a three-fold dis- 
crete spatial symmetry. Figure 12 shows a typical appli- 
cation of 1)14(1 , with the Newton descent search restricted 
to the configuration space. The initial loop. Fig. |2Ia), 
is a rather coarse initial guess. We fix arbitrarily the 
scaling A = 2.1, that is, we search for a cycle p of the 
fixed period Tp = 13.1947, with no constraint on the en- 
ergy. Figure l^fb) shows the cycle found by the Newton 
descent, with energy E ~ 0.1794, and the full discrete 
symmetry of the Hamiltonian. This cycle persists adia- 
batically for a small range of values of A; with A changed 
much, the Newton descent takes the same initial loop into 
other cycles. Figure EIc) verifies that the cost functional 
decreases exponentially with slope -2 throughout the 
r = [0, 10] integration interval, as predicted by ifT^ . The 
points get more and more sparse as r increases, because 
our numerical implementation adaptively chooses bigger 
and bigger step sizes 6t. 

In the Henon-Heiles case, the accelerations Qx , ay de- 
pend only on the configuration variables x, y. More gen- 
erally, the accelerations could also depend on x,y. Con- 
sider as an example the equations of motion for the re- 
stricted three-body problem [sif . 

, + X-1+ ^ 
X = 2y + .T-(l-^t) — 3 3 , 

y = -2i + 2/-(l-//)4-M4' (20) 

"^1 ^2 
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FIG. 2: The Henon-Heiles system in a chaotic region: (a) 
An initial loop L(0), and (b) the unstable periodic orbit p 
of period T = 13.1947 reached by the Newton descent H14|i . 
(c) The exponential decrease of the cost function, In(F^) ~ 
-2.0502 r + 6.0214. 



where t\ = ^(x^ fi)^ + y^, r2 = {x — 1 + /z)^ -I- y^. 
These equations describe the motion of a test particle in 
a rotating frame under the influence of the gravitational 
force of two heavy bodies with masses 1 and n <^ 1 fixed 
at {—fi,0) and (1 — /i,0) in the {x,y) coordinate frame. 
The stationary solutions of ((20|l are called the Lagrange 
points, corresponding to a circular motion of the test 
particle in phase with the rotation of the heavy bodies. 
The periodic solutions in the rotating frame correspond 
to periodic or quasi-periodic motion of the test particle 
in the inertial frame. Figure |31 shows an initial loop and 
the cycle to which it converges, in the rotating frame. 
Although the cycle looks simple, the Newton descent re- 
quires advancing in small St steps in order for the initial 
loop to converge to it. 

In order to successfully apply the Hamiltonian version 
of the Newton descent (|14|l , we have to ensure that the 
test particle keeps a finite distance from the origin. If 
a cycle passes very close to one of the heavy bodies. 
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FIG. 3: (a) An initial loop L{0), and (b) the unstable periodic 
orbit p of period Tp — 2.7365 reached by the Newton descent 
1141 1. the restricted three body problem 12UII in the chaotic 
regime, fx — 0.04. 

the acceleration can become so large that our scheme 
of uniformly distributing the loop points in time might 
fail to represent the loop faithfully. Another distribu- 
tion scheme is required in this case, for example, making 
the density of points proportional to the magnitude of 
acceleration. 



B. Periodic orbits of Kuramoto-Sivashinsky system 

The Kuramoto-Sivashinsky equation arises as an am- 
plitude equation for interfacial instability in a variety of 
contexts |32,|33- In 1-dimensional space, it reads 

Ut = ilL^)x - Uxx - VUxxxx, (21) 

where v is a, "super-viscosity" parameter which controls 
the rate of dissipation and {vF')x is the nonlinear con- 
vection term. As v decreases, the system undergoes a 
series of bifurcations, leading to increasingly turbulent, 
spatio-temporally chaotic dynamics. 

If we impose the periodic boundary condition u{t,x + 
2tt) = u{t,x) and choose to study only the odd solutions 
u(~x, t) = —u{x, t), u{x, t) can be expanded in a discrete 
spatial Fourier series [2J|, 

oo 

u{x,t)^i afcWe''", (22) 

k—~oo 



where a_fc = —a^ G K . In terms of the Fourier com- 
ponents, PDF (|21|) becomes an infinite ladder of ODEs, 

oo 

dk = (fc^ - vk'^)ak ~ X! "-niO-k-m ■ (23) 

rn— — oo 

In numerical simulations we work with the Galerkin trun- 
cations of the Fourier series since in the neighborhood of 
the strange attractor the magnitude of decreases very 
fast with k, high frequency modes playing a negligible 
role in the asymptotic dynamics. In this way Galerkin 
truncations reduce the dynamics to a finite but large 
number of ODEs. Wc work with d = 32 dimensions in 
our numerical calculations. In ref. p^ . multipoint shoot- 
ing has been successfully applied to obtain periodic orbits 
close to the onset of spatiotemporal chaos {v = 0.03). In 
this regime, our method is so stable that big time steps 
St can be employed even at the initial guesses, leading to 
extremely fast convergence. We attribute this robustness 
to the simplicity of the structure of the attractor at high 
viscosity values. 

The challenge comes with decreasing with the dy- 
namics turning more and more turbulent. Already at 
1/ = 0.015, the system is moderately turbulent and the 
phase space portraits of the dynamics reveal a complex 
labyrinth of "eddies" of different scales and orientations. 
While the highly unstable nature of orbits and intricate 
structure of the invariant set hinder applications of con- 
ventional cycle-search routines, in this setting our vari- 
ational method shines through. We design rather arbi- 
trary initial loops from numerical trajectory segments, 
and the calculation proceeds as before, except that now 
a small 6t has to be used initially to ensure numerical 
stability. Topologically different loops are very likely to 
result in different cycles, while some initial loop guesses 
my lead to local nonzero minima of the cost functional 
F^. As explained in Sect. II Vl in such cases the method 
diverges, and the search is restarted with a new inital 
loop guess. 

Two initial loop guesses are displayed in Fig. ^ along 
with the two periodic orbits detected by the Newton de- 
scent. In discretization of the initial loops, each point has 
to be specified in all d dimensions; here the coordinates 
{ai, 02} are picked so that topological similarity between 
initial and final loops is visually easy to identify. Other 
projections from d = 32 dimensions to subsets of 2 coor- 
dinates appear to make the identification harder, if not 
impossible. In both calculations, we molded segments 
of typical trajectories into smooth closed loops by the 
Fourier filtering method of Sect. lIVI As the desired orbit 
becomes longer and more complex, more sampling points 
are needed to represent the loop. We use N = 512 points 
to represent the loop in the (a)-(b) case and N = 1024 
points in the (c)-(d) case. The space-time evolution of 
u(x,t) for these two unstable spatio-temporally periodic 
solutions is displayed in Fig.|31 As u{x^ t) is antisymmet- 
ric on [— TT, tt], it suffices to display the solutions on the 
a; e [0, tt] interval. 
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FIG. 4: The Kuramoto-Sivashinsky system in a spatio- 
temporally turbulent regime (viscosity parameter v = 0.015, 
d = 32 Fourier modes truncation), (a) An initial guess Li, 
and (b) the periodic orbit pi of period T\ = 0.744892 reached 
by the Newton descent, (c) Another initial guess L2, and (d) 
the resulting periodic orbit p2 of period T2 — 1.184668. 



VI. DISCUSSION 

In order to cope with the difficulty of finding periodic 
orbits in high-dimensional chaotic flows, we have devised 
the Newton descent method, an infinitesimal step variant 
of the damped Newton-Raphson method. Our main re- 
sult is the PDE (|10|l which solves the variational problem 



0.6 



(a) 




(b) 

FIG. 5: Level plot of the space-time evolution u{x, t) for 
the two spatiotemporally periodic solutions of Fig. 2] (a) the 
evolution of pi, with start of a repeat after the cycle period 
Ti = 0.744892, and (b) one fuU period T2 = 1.184668 in the 
evolution of p2. 



of minimizing the cost functional This equation de- 
scribes the fictitious time t flow in the space of loops 
which decreases the cost functional at uniform exponen- 
tial rate (see ^VZ\). Variants of the method are presented 
for special classes of systems, such as Hamiltonian sys- 
tems. An efficient integration scheme for the PDE is de- 
vised and tested on the Kuramoto-Sivashinsky system, 
the Henon-Heiles system and the restricted three-body 
problem. 

Our method uses information from a large number of 
points in phase space, with the global topology of the 
desired cycle protected by insistence on smoothness and 
a uniform discretization of the loop. The method is quite 
robust in practice. 

The numerical results presented here are only a proof 
of principle. We do not know to what periodic orbit the 
flow H1U|) will evolve for a given dynamical system and a 
given initial loop. Empirically, the flow goes toward the 
"nearest" periodic orbit, with the largest topological re- 
semblance. Each particular application will still require 
much work in order to elucidate and enumerate relevant 
topological structures. The hope is that the short spatio- 
temporally periodic solutions revealed by the Newton 
descent searches will serve as the basic building blocks 
for systematic investigations of chaotic and perhaps even 
"turbulent" dynamics. 
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